---
title: "Gröbner bases are in 'qspray'"
author: "Stéphane Laurent"
date: '2023-07-14'
tags: R, maths
rbloggers: yes
output:
  md_document:
    variant: markdown
    preserve_yaml: true
  html_document:
    highlight: kate
    keep_md: no
highlighter: pandoc-solarized
---

There is something new in the **qspray** package: *Gröbner bases*.
Gröbner bases have many applications, and we will see one of them now.
My package **qspray** deals with multivariate polynomials with rational
coefficients.

Consider the following system of three polynomial equations:

`\begin{align}
x^2 + y + z & = 1 \\
x + y^2 + z & = 1 \\
x + y + z^2 & = 1
\end{align}`{=tex}

The goal is to find all solutions of this system, that is to say, the
simultaneous roots of the three polynomials $x^2 + y + z - 1$,
$x + y^2 + z - 1$, $x + y + z^2 - 1$. Not easy, is it?

*Fact:* a Gröbner basis of (the ideal generated by) a list of
polynomials is a list of polynomials with the same simultaneous roots.

``` r
library(qspray)
x <- qlone(1); y <- qlone(2); z <- qlone(3)
f1 <- x^2 + y + z - 1
f2 <- x + y^2 + z - 1
f3 <- x + y + z^2 - 1
gb <- groebner(list(f1, f2, f3)) # Gröbner basis of f1, f2, f3
( gbstrings <- lapply(gb, prettyQspray, vars = c("x", "y", "z")) )
## [[1]]
## [1] "y + x + z^2 - 1"
## 
## [[2]]
## [1] "((-9)*z^3)/2 - 2*z^5 - z^6 + 6*z^4 - y + z + z^7/2 + y^2"
## 
## [[3]]
## [1] "y*z^2 - z^2/2 + z^4/2"
## 
## [[4]]
## [1] "z^6 - z^2 + 4*z^3 - 4*z^4"
```

Do you see something nice? The fourth polynomial,
$z^6 - z^2 + 4z^3 - 4z^4$, is univariate: it depends only on the
variable $z$. We could numerically solve it with the R function
`polyroot`, but we want exact roots. It is easy to see that $0$ and $1$
are two roots of this polynomial. So we can divide this polynomial by
$z(z-1)$, and it will remain a factor of degree $4$. But I know, you are
lazy. Me too. That's why I prefer to use the **Ryacas** package to
perform this task:

``` r
library(Ryacas)
yac_str("Factor(z^6 - z^2 + 4*z^3 - 4*z^4)")
## [1] "(2*z+z^2-1)*(z-1)^2*z^2"
```

So it remains to find the roots of $2z + z^2 - 1$. We use **Ryacas**
again, as we're lazy:

``` r
yac_str("PSolve(2*z + z^2 - 1, z)")
## [1] "{(Sqrt(8)-2)/2,(-(Sqrt(8)+2))/2}"
```

There are two roots: $\sqrt{2} - 1$ and $-\sqrt{2} - 1$.

But let's consider the two friendly roots first: $0$ and $1$. If you
plug $z=0$ and $z=1$ in the three other polynomials of the Gröbner
basis, you easily find three simultaneous roots: $(1, 0, 0)$,
$(0, 1, 0)$ and $(0, 0, 1)$.

You should do this exercise, because now we will plug $z=\sqrt{2}-1$ and
$z=-\sqrt{2}-1$, and this will be less easy. In fact I will only treat
$z=\sqrt{2}-1$, and leave the other root as another exercise for you.

Let's go. We plug $z=\sqrt{2}-1$ in the third polynomial. Clearly, we
will obtain an univariate polynomial in $y$, and we will solve it. Or
rather, **Ryacas** will solve it:

``` r
yac_str(sprintf("p3 := %s", gbstrings[[3L]]))
## [1] "y*z^2-z^2/2+z^4/2"
yac_str("PSolve(Eliminate(z, Sqrt(2)-1, p3), y)")
## [1] "-((Sqrt(2)-1)^4/2-(Sqrt(2)-1)^2/2)/(Sqrt(2)-1)^2"
```

Yacas is not always nice. It does not simplify this expression. So let's
look at its numerical approximation:

``` r
yac_str("N(PSolve(Eliminate(z, Sqrt(2)-1, p3), y))")
## [1] "0.4142135623"
```

We recognize $\sqrt{2}-1$. Now we plug $z=\sqrt{2}-1$ into the second
polynomial and we solve the resulting polynomial of degree $2$ in $y$:

``` r
yac_str(sprintf("p2 := %s", gbstrings[[2L]]))
## [1] "((-9)*z^3)/2-2*z^5-z^6+6*z^4-y+z+z^7/2+y^2"
yac_str("N(PSolve(Eliminate(z, Sqrt(2)-1, p2), y))")
## [1] "{0.5857864376,0.4142135623}"
```

There are two roots: $y=\sqrt{2}-1$ and the other one has no importance,
since we didn't find it before.

So finally we just have to plug $y=\sqrt{2}-1$ and $z=\sqrt{2}-1$ into
the first polynomial and we solve the resulting univariate polynomial in
$x$:

``` r
yac_str(sprintf("p1 := %s", gbstrings[[1L]]))
## [1] "y+x+z^2-1"
yac_str(
  "N(PSolve(Eliminate(y, Sqrt(2)-1, Eliminate(z, Sqrt(2)-1, p1)), x))"
)
## [1] "0.4142135623"
```

Again? Well, we find $x = \sqrt{2}-1$. Finally we found one simultaneous
root: $$(x, y, z) = (\sqrt{2}-1, \sqrt{2}-1, \sqrt{2}-1).$$ Now I leave
you with the promised exercise: treat the case $z=-\sqrt{2}-1$. I give
you the solution, you will find:
$$(x, y, z) = (-\sqrt{2}-1, -\sqrt{2}-1, -\sqrt{2}-1).$$

Isn't it impressive? I started to learn about Gröbner bases only a few
days ago, and I still have to learn. I hope some applications will be
implemented in **qspray** in the future, and/or exposed on this blog.

The new version of **qspray** is not on CRAN yet. If you are in a hurry,
you can install the development version hosted [on
Github](https://github.com/stla/qspray).